% test spike-slab logistic model
addpath ../

ind = 1000;
npheno = 1;
nprobe = 10000;
causal = 100;

X = randn(ind,nprobe,'single');

effect_phe_true = randn(npheno,causal,'single');
pheno = single(1./(1+exp( - X(:,1:causal)*effect_phe_true' + 0.01*randn(ind,npheno,'single'))));

prior_var = mean(pheno(:).*(1-pheno(:)));

for pp = [-10, -8, -7, -5, -4, -3, -2, -1],

    prior = pp*ones(size(X,2),1,'single');
    [effect,effect_var,z,llik,elbo] = update_logit_spikeslab_model(pheno,X,prior,prior_var,100);
    fprintf('llik = %f, elbo = %f, Pr(Z) = %f\n',llik(end),elbo,mean(z));

end

figure(1); semilogx((effect./sqrt(effect_var))');
figure(3); semilogx(z);
figure(4); plot(llik)






